State Transition Tensors for Continuous-Thrust Control of Three-Body Relative Motion

Open AccessEngineering NotesState Transition Tensors for Continuous-Thrust Control of Three-Body Relative MotionJackson Kulik, William Clark and Dmitry SavranskyJackson KulikCornell University, Ithaca, New York 14850*Ph.D. Candidate, Center Applied Mathematics; .Search more papers by this author, ClarkCornell 14850†Visiting Assistant Professor, Mathematics.Search author SavranskyCornell 14850‡Associate Mechanical Aerospace Engineering.Search authorPublished Online:21 May 2023https://doi.org/10.2514/1.G007311SectionsRead Now ToolsAdd to favoritesDownload citationTrack citations ShareShare onFacebookTwitterLinked InRedditEmail AboutI. IntroductionTrajectory mission design problems often require bilevel optimizations determine optimal trajectories that balance goals with operational constraints. For example, in Refs. [1,2] authors deal doubly formation flight behavior the observation schedule maximizes science yield while minimizing control costs coordinated a telescope starshade Sun–Earth space. To meet these objectives, traveling-salesman-type problem is posed which stars observe on what schedule. A proposed itinerary judged partially basis fuel maneuver accordingly. The cost one leg journey may be assessed solving continuous-thrust, problem. Thus, computation multiple continuous-thrust just step evaluating objective function when considering overall scheduling Efficient algorithms solve relative context three-body dynamics are necessary facilitate solutions variety optimization problems. In particular, work will primarily focus not only computing control, but also quickly approximating control.The satellite another has been extensively studied two-body [3–9]. These methods generally rely analytical form state transition matrix Clohessy–Wiltshire or Tschauner–Hempel equations, available On other hand, impulsive formations well [10–13] some recent advances continuous control. Franzini et al. [14] examined applied motion elliptical using an adjoint method. While methodology improves fidelity speed computation, numerical integration still required each time new unsuitable applications high computational subproblem render overarching tractable. Park [15] presented method generating functions rendezvous problems, generalizes use where no known. approach allows solution arbitrary boundary conditions after single function, must recomputed from scratch whenever initial epoch adjusted. Similarly, Boone Mcmahon [16] achieved low-thrust differential dynamic programming tensors obtain feedback laws near reference orbit over specified flight. Their maneuvers allowed extended calculate flight, neighborhood two store temporal derivatives Taylor approximations [17,18]. However, do outside small region.We aim approximately epoch, without performing online. We seek limit storage requirements precomputed data, fitness metrics related integration. accomplish goals, we begin reviewing properties first- second-order variational equations. From there, formulate basic Finally, outline algorithm precompute interpolate data associated trajectory dynamical system states costates. Our interpolation scheme specially motivated power series equations employs variable order entry tensor being interpolated. interpolated linear costates few products twice size then allow approximation product involving tensor. Additionally, used Newton iteration find accurate than would information. Examples around Halo evaluate error approaches setting.A similar precomputation was taken previous space-telescope working conjunction [19]. Another interpolating uses regularization distribute points effectively Chebyshev polynomials perform [20]. This propose model Gateway mission, employed calculating maneuvers. Much present paper could performed within framework; however, explore leveraging cocycle resulting designed low memory coarse like emphasize before deriving our procedure intended analysis development ground tailored onboard satellite.II. Variational EquationsThe derived Ref. [21]. notion explored [22,23]. describe commonly known both component build intuition tensor, most clearly described form.Given autonomous Rn, vector x evolves according ordinary ddtx=F(x),x(0)=x0(1)The flow map defined such ddtφt(x)=F(φt(x)),φ0(x)=x(2)The possesses semigroup property φt○φs=φt+s(3)The Jacobian yields (STM) Φ(t,0) given 0 t. adopt indexing rows i columns j: Φji(t,0)=∂φti(x)∂xj(4)where upper index refers ith output. Exchanging spatial (assuming F derivatives) applying chain rule n2 first-order equations: dΦ(t,0)dt=∂F(x)∂xΦ(t,0),Φ(0,0)=In(5)where gives n-by-n identity matrix. Alternatively, components summation respect l understood δji Kronecker delta, dΦji(t,0)dt=∂Fi(x)∂xlΦjl(t,0),Φji(0,0)=δji(6)The recognizable stems application property: Φ(t2,t0)=Φ(t2,t1)Φ(t1,t0)(7)or, Φji(t2,t0)=Φli(t2,t1)Φjl(t1,t0)(8)Note relationship Eq. (8) commutative. Given knowledge STM entire interval stage, latter subinterval can calculated inverse Φ−1(t1,t0): Φ(t2,t1)=Φ(t2,t0)Φ−1(t1,t0)(9)Moving define (2,1)-state (STT) Ψ(t,0): Ψj,ki(t,0)=∂2φti(x)∂xj∂xk(10)Applying (6), n3 (11) depend values Ψ Φ: dΨj,ki(t,0)dt=∂2Fi(x)∂xl∂xqΦjl(t,0)Φkq(t,0)+∂Fi(x)∂xlΨj,kl(t,0),Ψj,ki(0,0)=0(11)The generalization comes differentiating (8): Ψj,ki(t2,t0)=Ψl,mi(t2,t1)Φjl(t1,t0)Φkm(t1,t0)+Φli(t2,t1)Ψj,kl(t1,t0)(12)Note made [24], contraction second copy left out. STTs along subinterval, STT finding (9): Ψj,ki(t2,t1)=[Ψl,mi(t2,t0)−Φqi(t2,t1)Ψl,mq(t1,t0)](Φ−1)jl(t1,t0)⁢(Φ−1)km(t1,t0)(13)This exact requires level approximate [20] map. With STT, perturbation truncated series: φti(x+δx)≈φti(x)+Φji(t,0)δxj+12Ψjki(t,0)δxjδxk(14)Variational describing sensitivity fixed parameters [25]. typically up first order, because they become increasingly complex as mixed partial appear at higher orders. avoid complications, note simply variables augmented have zero derivatives. sensitivities handled specific treating separately variables. potentially less efficient separately, additional multiplications additions conditions.III. Interpolating Short-Time State TensorsWe examine Φ(αt,0) α∈[0,1] Φ(t,0). straightforward way entrywise between time, Φ(0,0) (the matrix), original from, Φji(αt,0)≈In+α(Φji(t,0)−In)(15)Error (15) quadratic span Here, same number operations interpolation, greater accuracy. exploits structure underlying produce them. spans assuming stays constant trajectory: dΦ(t,0)dt≈∂F(x0)∂xΦ(t,0),Φ(0,0)=In(16)The equation exponential Φ(t,0)≈eAt=∑k=0∞tkk!Ak(17)where A=(∂F(x0)/∂x) A0=In convenience. Ak denotes kth A. approximated Φ(αt,0)≈eAαt=∑k=0∞(αt)kk!Ak(18)In components, tΦji(αt,0)≈In+αPji(Φji(t,0)−In)(19)where Pji=min{p∈N|(Ap)ji≠0}(20)That is, terms multiplied powers α determined nonzero term nonzero. energy circular restricted (CR3BP) later Eqs. (31) (39) consists mostly zeros. As itself repeatedly, fewer remain zero. Figure 1 depicts P CR3BP generic (corresponding undergoing natural motion). White regions graphical representation those undefined unnecessary.Next, consider Ψ(t,0) Ψ(αt,0). written unknown coefficient BN Ψj,ki(t,0)=∑N=0∞tNN!(BN)j,ki(21)For tΨj,ki(αt,0)≈αQj,kiΨj,ki(t,0)(22)where Qj,ki=min{p∈N|(Bp)j,ki≠0}(23)One BN; simplest Qj,ki numerically integrating short t t/2 rounding logarithm their quotient: Qj,ki=round(log2(Ψj,ki(t,0)Ψj,ki(t/2,0)))(24)We leading interpolation. Note expressions Φ(αt,0), nearly interval, Φ(t0+αt,t0)=Φ(αt,0) t0 if t0+αt≤t. true Ψ. primary benefit (19) multiplication per means take almost little manage higher-order entries. Only need online desired value Q computed system, regardless orbit.Fig. represents lowest corresponding nonzero.We functions. expect significantly outperform STT. compare short-time [Eq. (15)] (19)], refer Sec. VII comparison schemes example. errors orders magnitude lower instead summarized machinery algorithm. Next, under framework equations.IV. Optimal ProblemWe unconstrained thrust outlined [1,26]. applicable metric use, tractable treatment bounded possess discontinuous mass costate stemming switching bang-bang [26].The governing (1). paper, specifically analyze (39), though satisfy 12 final states, six-dimensional vectors consisting three-dimensional position velocity x=[rTvT]T: x(t0)=x0,x(tf)=xf(25)Simultaneously, minimize integral u: J=∫t0tf12uTu dt(26)The Hamiltonian (1) H=12uTu+λTF+pTu(27)where λ collection p=(λ4,λ5,λ6)T. u last three ddtλ=−(∂F(x)∂x)Tλ(28)A two-point arises (28) (25). Solution effort u=−p system. these, (26) integrated performance evaluated. Typically, repeated part shooting indirect generalized Legendre–Clebsch should checked local optimality: ∂2H∂u2>0(29)We attempt ahead needs solved By assume occurs trajectory.V. Approximate ControlAssume uncontrolled x(0)=x0, following until it reaches x(tf)=xf tf. chief set dynamics. interested controlling nearby deputy x0+δx0, so arrives xf+δxf δx(t) difference states.We now turn metric. Consider y=[xTλT]T,z=[xTλTJ]T(30)We I, z: ddtz=G(z)=[(F(x)T+[0uT])−λT(∂F(x)∂x)12uTu]T(31)In about applied. case, all (λ0=0). Φ optimally controlled vicinity trajectory. even stay trajectory, variations reflected since six G(z) linearly u, δx0 δxf manner: δλ0≈(Φλx(tf,t0))−1(δxf−Φxx(tf,t0)δx0)(32)where Φba(tf,t0)=∂a(tf)∂b(t0)(33)However, J(tf)=Jf cannot (it approximated, trivially zero). this, Jf≈12Ψy,yJ(t0,tf)δy0δy0(34)with Ψb,ca(tf,t0)=∂2a(tf)∂b(t0)∂c(t0)(35)Since J inherently quantity, ΦaJ=δaJ, δaJ delta. consequence fact (12)] independent tensor: Ψj,kJ(t2,t0)=Ψl,mJ(t2,t1)Φjl(t1,t0)Φkm(t1,t0)+Ψj,kJ(t1,t0)(36)Additionally, improve δλ0 (32) Newton’s (37) δxf: δxf≈Φyxδy0+12Ψy,yxδy0δy0(37)where Φλx derivative scheme. globally minimal x(0)=x0 satisfies any effort. condition satisfied further, (∂2H/∂u2)>0 sufficiently nonsingular Φλx(tf,t0), ensuring minimum (26).To summarize, t0, tf, compute (5) field 31. case 6-dimensional, 13-dimensional costates, cost. specified, directly iteration, respectively. solves problem, order. propagated (34). manner, close trajectory.VI. Precomputation DataTo expand orbit, develop data.We detail phase, phase solved. assumes range. periodic range assumed infinite. discretization 2m m∈N chosen, period broken into intervals equal length. Along integrated, Φ(Δm,j),Ψ(Δm,j) Δm,j stored Δm,j=(T0+Tf−T02mj,T0+Tf−T02m(j+1))(38)where T0 Tf give epochs consideration. zero, orbit.The (12) Φ(Δi,j),Ψ(Δi,j) decreasing m−1 0. m+1 different levels discretization, finest (or orbit) i=0. 2 multilevel orbit. it, Δi,j shown various indices half mth discretization. depicted finite horizons nonperiodic orbits modification. useful modeling higher-fidelity performs correction maneuvers, these. concludes algorithm, 2m+1(2n)2 floating point entries J, optionally 2m+1(2n)3 Ψyyy, n dimension 2n dimensions combined that, improved Ψy,yx (37). unlike Ψyyx(t2,t0) depends smaller Ψyyλ(t1,t0) (12). own sake. relatively fine (m=7 giving slightly single-day precision six-month L2 orbit), information double-precision format, amount 4 megabytes, fits manageably RAM modern desktop computer. Cutting does pertain usage [since (36) relies portions tensor] leads hundreds kilobytes. future, storing sparse directional rather whole cut accuracy [27]. Fourier possible reduce further here [28]. cubic spline suggested [20].For stage many succession. t0,tf Φ(tf,t0),Ψ(tf,t0) piecing together tensors.Fig. notional halfway labeled.Specifically, i* found there exists j* Δi*,j*⊆(t0,tf). here, process right disconnected inner removed, (t0,tf)∖Δi*,j*. running kept conditions. continues i=m, O(m) place Φ(tf′,t0′),Ψ(tf′,t0′), 0≤t0′−t0,tf−tf′≤(Tf−T0)2−m.We Φ(t0′,t0),Ψ(t0′,t0) Φ(tf,tf′),Ψ(tf,tf′) augment Φ(tf′,t0′),Ψ(tf′,t0′) come Φ(tf,t0),Ψ(tf,t0). (22), α=(t0′−t0)/(‖Δm,j‖) (t0,t0′)⊆Δm,j α=(tf−tf′)/(‖Δm,j‖) (tf′,tf)⊆Δm,j.With Φ(tf,t0) Ψ(tf,t0) obtained O(mn3) O(mn4) operations, respectively (due equations), dominated n-dimensional O(n3) complexity. Any future times O(n2) long relevant LU factorization (32). estimate (34) operations.A question how choose m. choice vary based tolerance error. tune m, smallest subintervals, compared against subinterval. norm normalized forward computations STM. increased acceptable level. check whether m absolute tolerances horizons, comparisons full calculations counterparts windows chosen